Chapter 4 Community composition

load("data/data.Rdata")

4.1 Taxonomy overview

4.1.1 Stacked barplot

# Merge data frames based on sample
merged_data<-genome_counts_filt %>%
  mutate_at(vars(-genome),~./sum(.)) %>% #apply TSS normalisation
  pivot_longer(-genome, names_to = "sample", values_to = "count") %>% #reduce to minimum number of columns
  left_join(., genome_metadata, by = join_by(genome == genome)) %>% #append genome metadata
  left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>% #append sample metadata
  filter(count > 0) #filter 0 counts

# Create an interaction variable for time_point and sample
merged_data$interaction_var <- interaction(merged_data$sample, merged_data$time_point)

ggplot(merged_data, aes(x=interaction_var,y=count, fill=phylum, group=phylum)) + #grouping enables keeping the same sorting of taxonomic units
    geom_bar(stat="identity", colour="white", linewidth=0.1) + #plot stacked bars with white borders
    scale_fill_manual(values=phylum_colors) +
    facet_nested(. ~ type,  scales="free") + #facet per day and treatment
    guides(fill = guide_legend(ncol = 1)) +
    theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1),
          axis.title.x = element_blank(),
          panel.background = element_blank(),
          panel.border = element_blank(),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          axis.line = element_line(linewidth = 0.5, linetype = "solid", colour = "black")) +
   labs(fill="Phylum",y = "Relative abundance",x="Samples")+
  scale_x_discrete(labels = function(x) gsub(".*\\.", "", x)) +
  facet_wrap(~type, scales = "free", labeller = as_labeller(function(label) gsub(".*\\.", "", label))) #only show time_point label

4.1.2 Phylum relative abundances

phylum_summary <- genome_counts_filt %>%
  mutate_at(vars(-genome),~./sum(.)) %>% #apply TSS normalisation
  pivot_longer(-genome, names_to = "sample", values_to = "count") %>%
  left_join(sample_metadata, by = join_by(sample == Tube_code)) %>%
  left_join(genome_metadata, by = join_by(genome == genome)) %>%
  group_by(sample,phylum) %>%
  summarise(relabun=sum(count))

phylum_summary %>%
    group_by(phylum) %>%
    summarise(mean=mean(relabun, na.rm=TRUE),sd=sd(relabun, na.rm=TRUE)) %>%
    arrange(-mean) %>%
    tt()
tinytable_lqmm6ff62lu6z187s5h8
phylum mean sd
p__Bacteroidota 0.380806816 0.202118047
p__Bacillota_A 0.253813952 0.163904042
p__Bacillota 0.118079729 0.146679823
p__Pseudomonadota 0.096901291 0.160581847
p__Campylobacterota 0.053050348 0.092998130
p__Verrucomicrobiota 0.029867960 0.073128363
p__Desulfobacterota 0.023580475 0.036482484
p__Chlamydiota 0.010572704 0.059690269
p__Fusobacteriota 0.010482132 0.028311738
p__Cyanobacteriota 0.009206465 0.016492133
p__Bacillota_C 0.004713164 0.006645703
p__Spirochaetota 0.004009680 0.012308028
p__Bacillota_B 0.002465465 0.004927667
p__Actinomycetota 0.001235741 0.006346852
p__Elusimicrobiota 0.001214079 0.006501752
phylum_arrange <- phylum_summary %>%
    group_by(phylum) %>%
    summarise(mean=mean(relabun)) %>%
    arrange(-mean) %>%
    select(phylum) %>%
    pull()

phylum_summary %>%
    filter(phylum %in% phylum_arrange) %>%
    mutate(phylum=factor(phylum,levels=rev(phylum_arrange))) %>%
    ggplot(aes(x=relabun, y=phylum, group=phylum, color=phylum)) +
        scale_color_manual(values=phylum_colors[rev(phylum_arrange)]) +
        geom_jitter(alpha=0.5) + 
        theme_minimal() + 
        theme(legend.position="none") +
        labs(y="Phylum",x="Relative abundance")

4.2 Taxonomy boxplot

4.2.1 Family

family_summary <- genome_counts_filt %>%
  mutate_at(vars(-genome),~./sum(.)) %>% #apply TSS normalisation
  pivot_longer(-genome, names_to = "sample", values_to = "count") %>% #reduce to minimum number of columns
  left_join(sample_metadata, by = join_by(sample == Tube_code)) %>% #append sample metadata
  left_join(., genome_metadata, by = join_by(genome == genome)) %>% #append genome metadata
  group_by(sample,family) %>%
  summarise(relabun=sum(count))

family_summary %>%
    group_by(family) %>%
    summarise(mean=mean(relabun, na.rm=TRUE),sd=sd(relabun, na.rm=TRUE)) %>%
    arrange(-mean) %>%
    tt()
tinytable_9i460mg5yncmsnzhsdue
family mean sd
f__Bacteroidaceae 2.215146e-01 0.1392048812
f__Lachnospiraceae 1.406192e-01 0.1085410432
f__Tannerellaceae 1.028077e-01 0.0799387840
f__Helicobacteraceae 5.258990e-02 0.0925721039
f__Mycoplasmoidaceae 3.694491e-02 0.0756423660
f__Erysipelotrichaceae 3.502391e-02 0.0451216980
f__UBA3700 3.418595e-02 0.0557485850
f__Marinifilaceae 2.774048e-02 0.0271999075
f__ 2.772765e-02 0.0838376915
f__DTU072 2.700317e-02 0.0529029477
f__Rikenellaceae 2.696900e-02 0.0464477971
f__Enterobacteriaceae 2.691677e-02 0.0918386902
f__Coprobacillaceae 2.551773e-02 0.0892179821
f__Desulfovibrionaceae 2.358047e-02 0.0364824837
f__Ruminococcaceae 1.859938e-02 0.0429483670
f__LL51 1.759437e-02 0.0687507595
f__Rhizobiaceae 1.596807e-02 0.0770633474
f__UBA3830 1.578817e-02 0.0454016753
f__Akkermansiaceae 1.227359e-02 0.0316838689
f__Chlamydiaceae 1.057270e-02 0.0596902690
f__Fusobacteriaceae 1.048213e-02 0.0283117384
f__CAG-239 9.141331e-03 0.0150992357
f__Enterococcaceae 8.420523e-03 0.0466580889
f__Gastranaerophilaceae 7.728770e-03 0.0144404459
f__Oscillospiraceae 6.643218e-03 0.0078105983
f__UBA1997 6.378408e-03 0.0309668631
f__Streptococcaceae 6.366441e-03 0.0342174908
f__UBA1242 4.673359e-03 0.0153730061
f__Brevinemataceae 4.009680e-03 0.0123080282
f__Acutalibacteraceae 3.374550e-03 0.0109560112
f__UBA660 3.196511e-03 0.0139558713
f__Clostridiaceae 2.842205e-03 0.0171339212
f__RUG11792 2.817729e-03 0.0251127757
f__Peptococcaceae 2.465465e-03 0.0049276669
f__MGBC116941 2.160169e-03 0.0093827415
f__Acidaminococcaceae 1.943804e-03 0.0050321931
f__CAG-508 1.807978e-03 0.0064175379
f__Anaerovoracaceae 1.569531e-03 0.0036471440
f__Moraxellaceae 1.485415e-03 0.0097446552
f__RUG14156 1.477695e-03 0.0045847831
f__Staphylococcaceae 1.361709e-03 0.0050863070
f__Elusimicrobiaceae 1.214079e-03 0.0065017516
f__CAG-288 9.490865e-04 0.0060198775
f__Anaerotignaceae 8.989745e-04 0.0040469504
f__CALVMC01 7.516846e-04 0.0043609744
f__Eggerthellaceae 6.407882e-04 0.0021266467
f__Massilibacillaceae 6.235073e-04 0.0016350480
f__Mycobacteriaceae 5.949530e-04 0.0060400054
f__UBA1820 4.736030e-04 0.0013057353
f__Arcobacteraceae 4.604457e-04 0.0050324221
f__CAG-274 4.519746e-04 0.0022028477
f__Burkholderiaceae_C 3.699430e-04 0.0048092594
f__Muribaculaceae 3.617894e-04 0.0009776409
f__UBA932 3.419549e-04 0.0011583178
f__Hepatoplasmataceae 2.989107e-04 0.0038858387
f__Rhodobacteraceae 2.959092e-04 0.0038468195
f__Weeksellaceae 2.771627e-04 0.0031476408
f__Eubacteriaceae 1.646822e-04 0.0006729067
f__Sphingobacteriaceae 1.505775e-04 0.0012460017
f__Devosiaceae 1.489995e-04 0.0015094318
f__Pumilibacteraceae 1.277418e-04 0.0007646754
f__WRAU01 9.603359e-05 0.0012484367
f__Peptostreptococcaceae 2.287338e-05 0.0002973539
family_arrange <- family_summary %>%
    group_by(family) %>%
    summarise(mean=sum(relabun)) %>%
    arrange(-mean) %>%
    select(family) %>%
    pull()

family_summary %>%
    left_join(genome_metadata %>% select(family,phylum) %>% unique(),by=join_by(family==family)) %>%
    left_join(sample_metadata,by=join_by(sample==Tube_code)) %>%
    filter(family %in% family_arrange[1:20]) %>%
    mutate(family=factor(family,levels=rev(family_arrange[1:20]))) %>%
    filter(relabun > 0) %>%
    ggplot(aes(x=relabun, y=family, group=family, color=phylum)) +
        scale_color_manual(values=phylum_colors[-8]) +
        geom_jitter(alpha=0.5) + 
        facet_grid(.~type)+
        theme_minimal() + 
        labs(y="Family", x="Relative abundance", color="Phylum")

4.2.2 Genus

genus_summary <- genome_counts_filt %>%
  mutate_at(vars(-genome),~./sum(.)) %>% #apply TSS nornalisation
  pivot_longer(-genome, names_to = "sample", values_to = "count") %>% #reduce to minimum number of columns
  left_join(sample_metadata, by = join_by(sample == Tube_code)) %>% #append sample metadata
  left_join(genome_metadata, by = join_by(genome == genome)) %>% #append genome metadata
  group_by(sample,genus) %>%
  summarise(relabun=sum(count)) %>%
  filter(genus != "g__")

genus_summary %>%
    group_by(genus) %>%
    summarise(mean=mean(relabun, na.rm=TRUE),sd=sd(relabun, na.rm=TRUE)) %>%
    arrange(-mean) %>%
    tt()
tinytable_k1fvanpcuw1fz73njhys
genus mean sd
g__Bacteroides 1.352836e-01 0.0923254347
g__Parabacteroides 9.681305e-02 0.0801790006
g__Phocaeicola 6.935126e-02 0.0794715561
g__Mycoplasmoides 3.051289e-02 0.0753256189
g__Helicobacter_J 3.009004e-02 0.0594808174
g__Odoribacter 2.552307e-02 0.0267570773
g__Roseburia 2.384718e-02 0.0559439989
g__NHYM01 2.249986e-02 0.0796951054
g__Alistipes 2.208351e-02 0.0283591693
g__Coprobacillus 2.013859e-02 0.0878851683
g__Agrobacterium 1.596807e-02 0.0770633474
g__Akkermansia 1.227359e-02 0.0316838689
g__Fusobacterium_A 1.038903e-02 0.0283170720
g__Kineothrix 8.801109e-03 0.0409037874
g__Proteus 8.657984e-03 0.0681831315
g__Dielma 8.578108e-03 0.0090836958
g__CAG-95 8.109350e-03 0.0205019603
g__JAAYNV01 7.994602e-03 0.0195404994
g__UBA866 7.260325e-03 0.0293499698
g__Desulfovibrio 7.018155e-03 0.0211479806
g__Enterococcus 7.001300e-03 0.0456600395
g__Ureaplasma 6.432013e-03 0.0138012183
g__Lactococcus 6.366441e-03 0.0342174908
g__Lacrimispora 6.064761e-03 0.0098190433
g__Parabacteroides_B 5.994680e-03 0.0100174097
g__CALXRO01 5.765728e-03 0.0308524398
g__Citrobacter 5.717035e-03 0.0334547978
g__NSJ-61 5.560111e-03 0.0198839325
g__Breznakia 5.504528e-03 0.0237016548
g__Clostridium_AQ 5.326190e-03 0.0121695008
g__Salmonella 5.133583e-03 0.0184556360
g__Bilophila 4.983840e-03 0.0088456017
g__Hungatella_A 4.811802e-03 0.0095549574
g__MGBC136627 4.364796e-03 0.0163003355
g__Escherichia 4.188365e-03 0.0266100578
g__UMGS1251 4.159842e-03 0.0072716746
g__Hungatella 4.116388e-03 0.0190788279
g__Clostridium_Q 4.011996e-03 0.0052127439
g__Brevinema 4.009680e-03 0.0123080282
g__Thomasclavelia 3.902580e-03 0.0109042017
g__Scatousia 3.649941e-03 0.0102727788
g__Enterocloster 3.620384e-03 0.0047339282
g__Mailhella 3.612079e-03 0.0102470769
g__Copromonas 3.599368e-03 0.0050180473
g__Ventrimonas 3.513355e-03 0.0071233140
g__Caccovivens 3.343316e-03 0.0122194941
g__Lawsonia 3.289015e-03 0.0117180962
g__Fournierella 3.217226e-03 0.0062309921
g__Limenecus 3.163279e-03 0.0065731820
g__MGBC133411 3.042418e-03 0.0074576633
g__Mucinivorans 2.900095e-03 0.0373193950
g__Sarcina 2.842205e-03 0.0171339212
g__Acetatifactor 2.738138e-03 0.0058449483
g__Eisenbergiella 2.697412e-03 0.0068331645
g__Bacteroides_G 2.682722e-03 0.0346150378
g__CAJLXD01 2.633818e-03 0.0087495783
g__Blautia 2.558708e-03 0.0061407647
g__C-19 2.275515e-03 0.0048691968
g__Butyricimonas 2.217402e-03 0.0050081070
g__Velocimicrobium 2.201224e-03 0.0066764022
g__CAZU01 2.196385e-03 0.0065944622
g__MGBC116941 2.160169e-03 0.0093827415
g__Intestinimonas 2.081898e-03 0.0039382613
g__Negativibacillus 2.069077e-03 0.0055137501
g__Rikenella 1.985389e-03 0.0037067264
g__Phascolarctobacterium 1.943804e-03 0.0050321931
g__RGIG6463 1.789843e-03 0.0039682804
g__JALFVM01 1.742360e-03 0.0038623852
g__Oscillibacter 1.510459e-03 0.0024992661
g__Pseudoflavonifractor 1.502976e-03 0.0027706264
g__Acinetobacter 1.485415e-03 0.0097446552
g__Citrobacter_A 1.392923e-03 0.0060348874
g__Staphylococcus 1.361709e-03 0.0050863070
g__RGIG4733 1.303790e-03 0.0040480788
g__UBA1436 1.214079e-03 0.0065017516
g__Lachnotalea 1.208197e-03 0.0049169585
g__Ruthenibacterium 1.202021e-03 0.0053632179
g__14-2 1.184643e-03 0.0096299356
g__Beduini 1.173950e-03 0.0025142208
g__Scatocola 1.121193e-03 0.0044975564
g__Faecisoma 1.085585e-03 0.0055596648
g__Enterococcus_A 1.083991e-03 0.0099150522
g__Faecimonas 9.880526e-04 0.0083079244
g__RGIG9287 9.638826e-04 0.0093228400
g__CAG-345 9.490865e-04 0.0060198775
g__Blautia_A 9.208028e-04 0.0029082304
g__RGIG1896 8.343610e-04 0.0062772562
g__CAG-269 7.982448e-04 0.0047176841
g__Marseille-P3106 7.938268e-04 0.0017618833
g__WRHT01 6.429563e-04 0.0026979819
g__Eggerthella 6.407882e-04 0.0021266467
g__IOR16 6.398942e-04 0.0020651222
g__Anaerotruncus 6.265445e-04 0.0016948517
g__RUG14156 6.151780e-04 0.0022147136
g__CHH4-2 6.145042e-04 0.0019997615
g__Corynebacterium 5.949530e-04 0.0060400054
g__Serratia_A 5.860616e-04 0.0076188009
g__CAG-56 4.915432e-04 0.0016341973
g__Merdimorpha 4.736030e-04 0.0013057353
g__MGBC140009 4.679334e-04 0.0024067320
g__CALURL01 4.635287e-04 0.0016737486
g__Aliarcobacter 4.604457e-04 0.0050324221
g__Scatenecus 4.520215e-04 0.0019760876
g__RGIG8482 4.399064e-04 0.0029753981
g__Enterobacter 4.073437e-04 0.0041317730
g__Klebsiella 4.054439e-04 0.0048910857
g__Caccenecus 3.941198e-04 0.0017802371
g__Alcaligenes 3.699430e-04 0.0048092594
g__Plesiomonas 3.633249e-04 0.0027105056
g__HGM05232 3.617894e-04 0.0009776409
g__JAHHSE01 3.592120e-04 0.0014900895
g__Egerieousia 3.419549e-04 0.0011583178
g__Emergencia 3.413630e-04 0.0017450341
g__Enterococcus_B 3.352316e-04 0.0022266910
g__Stoquefichus 3.026072e-04 0.0020503969
g__Hepatoplasma 2.989107e-04 0.0038858387
g__Paracoccus 2.959092e-04 0.0038468195
g__Moheibacter 2.771627e-04 0.0031476408
g__Scatomorpha 2.641015e-04 0.0010184339
g__UBA7185 2.434328e-04 0.0014558191
g__Eubacterium 1.646822e-04 0.0006729067
g__Sphingobacterium 1.505775e-04 0.0012460017
g__Devosia 1.489995e-04 0.0015094318
g__Anaerosporobacter 1.454112e-04 0.0012747262
g__Caccomorpha 1.383122e-04 0.0010540604
g__UBA2658 1.307451e-04 0.0007204965
g__Protoclostridium 1.277418e-04 0.0007646754
g__Angelakisella 1.268479e-04 0.0009221501
g__Cetobacterium_A 9.310217e-05 0.0008718575
g__Rahnella 6.470705e-05 0.0008411917
g__Peptostreptococcus 2.287338e-05 0.0002973539
genus_arrange <- genus_summary %>%
    group_by(genus) %>%
    summarise(mean=sum(relabun)) %>%
    filter(genus != "g__")%>%
    arrange(-mean) %>%
    select(genus) %>%
    mutate(genus= sub("^g__", "", genus)) %>%
    pull()

genus_summary %>%
    left_join(genome_metadata %>% select(genus,phylum) %>% unique(),by=join_by(genus==genus)) %>%
    left_join(sample_metadata,by=join_by(sample==Tube_code)) %>%
    mutate(genus= sub("^g__", "", genus)) %>%
    filter(genus %in% genus_arrange[1:20]) %>%
    mutate(genus=factor(genus,levels=rev(genus_arrange[1:20]))) %>%
    filter(relabun > 0) %>%
    ggplot(aes(x=relabun, y=genus, group=genus, color=phylum)) +
        scale_color_manual(values=phylum_colors[-c(3,4,6,8)]) +
        geom_jitter(alpha=0.5) + 
        facet_grid(.~type)+
        theme_minimal() + 
        labs(y="Family", x="Relative abundance", color="Phylum")

4.3 Alpha diversity

# Calculate Hill numbers
richness <- genome_counts_filt %>%
  column_to_rownames(var = "genome") %>%
  dplyr::select(where(~ !all(. == 0))) %>%
  hilldiv(., q = 0) %>%
  t() %>%
  as.data.frame() %>%
  dplyr::rename(richness = 1) %>%
  rownames_to_column(var = "sample")

neutral <- genome_counts_filt %>%
  column_to_rownames(var = "genome") %>%
  dplyr::select(where(~ !all(. == 0))) %>%
  hilldiv(., q = 1) %>%
  t() %>%
  as.data.frame() %>%
  dplyr::rename(neutral = 1) %>%
  rownames_to_column(var = "sample")

phylogenetic <- genome_counts_filt %>%
  column_to_rownames(var = "genome") %>%
  dplyr::select(where(~ !all(. == 0))) %>%
  hilldiv(., q = 1, tree = genome_tree) %>%
  t() %>%
  as.data.frame() %>%
  dplyr::rename(phylogenetic = 1) %>%
  rownames_to_column(var = "sample")

# Aggregate basal GIFT into elements
dist <- genome_gifts %>%
  to.elements(., GIFT_db) %>%
  traits2dist(., method = "gower")

functional <- genome_counts_filt %>%
  column_to_rownames(var = "genome") %>%
  dplyr::select(where(~ !all(. == 0))) %>%
  hilldiv(., q = 1, dist = dist) %>%
  t() %>%
  as.data.frame() %>%
  dplyr::rename(functional = 1) %>%
  rownames_to_column(var = "sample") %>%
  mutate(functional = if_else(is.nan(functional), 1, functional))

# Merge all metrics
alpha_div <- richness %>%
  full_join(neutral, by = join_by(sample == sample)) %>%
  full_join(phylogenetic, by = join_by(sample == sample)) %>%
  full_join(functional, by = join_by(sample == sample))

4.3.1 Wild samples

alpha_div %>%
  pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
  left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
  filter(time_point=="0_Wild") %>%
  mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic","functional"))) %>%
      ggplot(aes(y = value, x = species, group=species, color=species, fill=species)) +
      geom_boxplot(outlier.shape = NA) +
      geom_jitter(alpha=0.5) +
      scale_color_manual(name="species",
          breaks=c("Podarcis_liolepis","Podarcis_muralis"),
          labels=c("Podarcis liolepis","Podarcis muralis"),
          values=c("#e5bd5b", "#6b7398")) +
      scale_fill_manual(name="species",
          breaks=c("Podarcis_liolepis","Podarcis_muralis"),
          labels=c("Podarcis liolepis","Podarcis muralis"),
          values=c("#e5bd5b50", "#6b739850")) +
      facet_wrap(. ~ metric ) +
      coord_cartesian(xlim = c(1, NA)) +
      theme_classic() +
        theme(
    strip.background = element_blank(),
    panel.grid.minor.x = element_line(size = .1, color = "grey"),
    axis.title.x = element_blank(),
    axis.title.y = element_blank(),
    axis.text.x = element_text(angle = 45, hjust = 1),
    # Increase plot size
    plot.title = element_text(size = 10),
    axis.text = element_text(size = 8),
    axis.title = element_text(size = 8)
      )

4.3.2 Acclimation samples

alpha_div %>%
  pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
  left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
  filter(time_point=="1_Acclimation") %>%
  mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic","functional"))) %>%
      ggplot(aes(y = value, x = species, group=species, color=species, fill=species)) +
      geom_boxplot(outlier.shape = NA) +
      geom_jitter(alpha=0.5) +
      scale_color_manual(name="species",
          breaks=c("Podarcis_liolepis","Podarcis_muralis"),
          labels=c("Podarcis liolepis","Podarcis muralis"),
          values=c("#e5bd5b", "#6b7398")) +
      scale_fill_manual(name="species",
          breaks=c("Podarcis_liolepis","Podarcis_muralis"),
          labels=c("Podarcis liolepis","Podarcis muralis"),
          values=c("#e5bd5b50", "#6b739850")) +
      facet_wrap(. ~ metric) +
      coord_cartesian(xlim = c(1, NA)) +
      theme_classic() +
        theme(
    strip.background = element_blank(),
    panel.grid.minor.x = element_line(size = .1, color = "grey"),
    axis.title.x = element_blank(),
    axis.title.y = element_blank(),
    axis.text.x = element_text(angle = 45, hjust = 1),
    # Increase plot size
    plot.title = element_text(size = 10),
    axis.text = element_text(size = 8),
    axis.title = element_text(size = 8)
      )

4.3.3 Antibiotics samples

alpha_div %>%
  pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
  left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
  filter(time_point=="2_Antibiotics") %>%
  mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic","functional"))) %>%
      ggplot(aes(y = value, x = species, group=species, color=species, fill=species)) +
      geom_boxplot(outlier.shape = NA) +
      geom_jitter(alpha=0.5) +
      scale_color_manual(name="species",
          breaks=c("Podarcis_liolepis","Podarcis_muralis"),
          labels=c("Podarcis liolepis","Podarcis muralis"),
          values=c("#e5bd5b", "#6b7398")) +
      scale_fill_manual(name="species",
          breaks=c("Podarcis_liolepis","Podarcis_muralis"),
          labels=c("Podarcis liolepis","Podarcis muralis"),
          values=c("#e5bd5b50", "#6b739850")) +
      facet_wrap(. ~ metric) +
      coord_cartesian(xlim = c(1, NA)) +
      theme_classic() +
        theme(
    strip.background = element_blank(),
    panel.grid.minor.x = element_line(size = .1, color = "grey"),
    axis.title.x = element_blank(),
    axis.title.y = element_blank(),
    axis.text.x = element_text(angle = 45, hjust = 1),
    # Increase plot size
    plot.title = element_text(size = 10),
    axis.text = element_text(size = 8),
    axis.title = element_text(size = 8)
      )

4.3.4 Transplant_1 samples

alpha_div %>%
  pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
  left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
  filter(time_point=="3_Transplant1") %>%
  mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic","functional"))) %>%
      ggplot(aes(y = value, x = type, group=type, color=type, fill=type)) +
      geom_boxplot(outlier.shape = NA) +
      geom_jitter(alpha=0.5) +
      scale_color_manual(name="type",
          breaks=c("Control","Hot_control", "Treatment"),
          labels=c("Control","Hot control", "Treatment"),
          values=c("#76b183","#d57d2c", "#4477AA")) +
      scale_fill_manual(name="type",
          breaks=c("Control","Hot_control", "Treatment"),
          labels=c("Control","Hot control", "Treatment"),
          values=c("#76b18350","#d57d2c50", "#4477AA50")) +
      facet_wrap(. ~ metric) +
      coord_cartesian(xlim = c(1, NA)) +
      theme_classic() +
        theme(
    strip.background = element_blank(),
    panel.grid.minor.x = element_line(size = .1, color = "grey"),
    axis.title.x = element_blank(),
    axis.title.y = element_blank(),
    axis.text.x = element_text(angle = 45, hjust = 1),
    # Increase plot size
    plot.title = element_text(size = 10),
    axis.text = element_text(size = 8),
    axis.title = element_text(size = 8)
      )

4.3.5 Transplant_2 samples

alpha_div %>%
  pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
  left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
  filter(time_point=="4_Transplant2") %>%
  mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic","functional"))) %>%
      ggplot(aes(y = value, x = type, group=type, color=type, fill=type)) +
      geom_boxplot(outlier.shape = NA) +
      geom_jitter(alpha=0.5) +
      scale_color_manual(name="type",
          breaks=c("Control","Hot_control", "Treatment"),
          labels=c("Control","Hot control", "Treatment"),
          values=c("#76b183","#d57d2c", "#4477AA")) +
      scale_fill_manual(name="type",
          breaks=c("Control","Hot_control", "Treatment"),
          labels=c("Control","Hot control", "Treatment"),
          values=c("#76b18350","#d57d2c50", "#4477AA50")) +
      facet_wrap(. ~ metric) +
      coord_cartesian(xlim = c(1, NA)) +
      theme_classic() +
        theme(
    strip.background = element_blank(),
    panel.grid.minor.x = element_line(size = .1, color = "grey"),
    axis.title.x = element_blank(),
    axis.title.y = element_blank(),
    axis.text.x = element_text(angle = 45, hjust = 1),
    # Increase plot size
    plot.title = element_text(size = 10),
    axis.text = element_text(size = 8),
    axis.title = element_text(size = 8)
      )

4.3.6 Post-Transplant_1 samples

alpha_div %>%
  pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
  left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
  filter(time_point=="5_Post-FMT1") %>%
  mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic","functional"))) %>%
      ggplot(aes(y = value, x = type, group=type, color=type, fill=type)) +
      geom_boxplot(outlier.shape = NA) +
      geom_jitter(alpha=0.5) +
      scale_color_manual(name="type",
          breaks=c("Control","Hot_control", "Treatment"),
          labels=c("Control","Hot control", "Treatment"),
          values=c("#76b183","#d57d2c", "#4477AA")) +
      scale_fill_manual(name="type",
          breaks=c("Control","Hot_control", "Treatment"),
          labels=c("Control","Hot control", "Treatment"),
          values=c("#76b18350","#d57d2c50", "#4477AA50")) +
      facet_wrap(. ~ metric) +
      coord_cartesian(xlim = c(1, NA)) +
      theme_classic() +
        theme(
    strip.background = element_blank(),
    panel.grid.minor.x = element_line(size = .1, color = "grey"),
    axis.title.x = element_blank(),
    axis.title.y = element_blank(),
    axis.text.x = element_text(angle = 45, hjust = 1),
    # Increase plot size
    plot.title = element_text(size = 10),
    axis.text = element_text(size = 8),
    axis.title = element_text(size = 8)
      )

4.3.7 Post-Transplant_2 samples

alpha_div %>%
  pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
  left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
  filter(time_point=="6_Post-FMT2") %>%
  mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic","functional"))) %>%
      ggplot(aes(y = value, x = type, group=type, color=type, fill=type)) +
      geom_boxplot(outlier.shape = NA) +
      geom_jitter(alpha=0.5) +
      scale_color_manual(name="type",
          breaks=c("Control","Hot_control", "Treatment"),
          labels=c("Control","Hot control", "Treatment"),
          values=c("#76b183","#d57d2c", "#4477AA")) +
      scale_fill_manual(name="type",
          breaks=c("Control","Hot_control", "Treatment"),
          labels=c("Control","Hot control", "Treatment"),
          values=c("#76b18350","#d57d2c50", "#4477AA50")) +
      facet_wrap(. ~ metric) +
      coord_cartesian(xlim = c(1, NA)) +
      theme_classic() +
        theme(
    strip.background = element_blank(),
    panel.grid.minor.x = element_line(size = .1, color = "grey"),
    axis.title.x = element_blank(),
    axis.title.y = element_blank(),
    axis.text.x = element_text(angle = 45, hjust = 1),
    # Increase plot size
    plot.title = element_text(size = 10),
    axis.text = element_text(size = 8),
    axis.title = element_text(size = 8)
      )

4.4 Beta diversity

beta_q0n <- genome_counts_filt %>%
  column_to_rownames(., "genome") %>%
  hillpair(., q = 0)

beta_q1n <- genome_counts_filt %>%
  column_to_rownames(., "genome") %>%
  hillpair(., q = 1)

beta_q1p <- genome_counts_filt %>%
  column_to_rownames(., "genome") %>%
  hillpair(., q = 1, tree = genome_tree)

beta_q1f <- genome_counts_filt %>%
  column_to_rownames(., "genome") %>%
  hillpair(., q = 1, dist = dist)

4.5 Permanovas

4.5.0.1 Load required data

meta <- column_to_rownames(sample_metadata, "Tube_code")

4.5.1 1. Are the wild populations similar?

4.5.1.1 Wild: P.muralis vs P.liolepis

wild <- meta %>%
  filter(time_point == "0_Wild")

# Create a temporary modified version of genome_counts_filt
temp_genome_counts <- transform(genome_counts_filt, row.names = genome_counts_filt$genome)
temp_genome_counts$genome <- NULL

wild.counts <- temp_genome_counts[, which(colnames(temp_genome_counts) %in% rownames(wild))]
identical(sort(colnames(wild.counts)), sort(as.character(rownames(wild))))

wild_nmds <- sample_metadata %>%
  filter(time_point == "0_Wild")

4.5.1.2 Number of samples used

[1] 28
beta_div_richness_wild<-hillpair(data=wild.counts, q=0)
beta_div_neutral_wild<-hillpair(data=wild.counts, q=1)
beta_div_phylo_wild<-hillpair(data=wild.counts, q=1, tree=genome_tree)
beta_div_func_wild<-hillpair(data=wild.counts, q=1, dist=dist)

4.5.1.3 Richness

Df SumOfSqs R2 F Pr(>F)
species 1 1.631953 0.207198 6.795072 0.001
Residual 26 6.244347 0.792802 NA NA
Total 27 7.876300 1.000000 NA NA

4.5.1.4 Neutral

Df SumOfSqs R2 F Pr(>F)
species 1 2.018797 0.2566342 8.976049 0.001
Residual 26 5.847641 0.7433658 NA NA
Total 27 7.866438 1.0000000 NA NA

4.5.1.5 Phylogenetic

Df SumOfSqs R2 F Pr(>F)
species 1 0.3786052 0.2108638 6.947419 0.001
Residual 26 1.4168908 0.7891362 NA NA
Total 27 1.7954960 1.0000000 NA NA

4.5.1.6 Functional

Df SumOfSqs R2 F Pr(>F)
species 1 0.0787916 0.1594096 4.930642 0.068
Residual 26 0.4154800 0.8405904 NA NA
Total 27 0.4942716 1.0000000 NA NA
beta_q0n_nmds_wild <- beta_div_richness_wild$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE, trace=FALSE) %>%
                vegan::scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(wild_nmds, by = join_by(sample == Tube_code))

beta_q1n_nmds_wild <- beta_div_neutral_wild$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE, trace=FALSE) %>%
                vegan::scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(wild_nmds, by = join_by(sample == Tube_code))

beta_q1p_nmds_wild <- beta_div_phylo_wild$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                vegan::scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(wild_nmds, by = join_by(sample == Tube_code))

beta_q1f_nmds_wild <- beta_div_func_wild$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                vegan::scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(wild_nmds, by = join_by(sample == Tube_code))

4.5.2 2. Do the antibiotics work?

4.5.2.1 Acclimation vs antibiotics

treat <- meta  %>% #antibiotics samples giving error when plotting
  filter(time_point == "1_Acclimation" | time_point == "2_Antibiotics")

temp_genome_counts <- transform(genome_counts_filt, row.names = genome_counts_filt$genome)
temp_genome_counts$genome <- NULL

treat.counts <- temp_genome_counts[,which(colnames(temp_genome_counts) %in% rownames(treat))]
identical(sort(colnames(treat.counts)),sort(as.character(rownames(treat))))

treat_nmds <- sample_metadata %>%
  filter(time_point == "1_Acclimation" | time_point == "2_Antibiotics")

4.5.2.2 Number of samples used

[1] 53
beta_div_richness_treat<-hillpair(data=treat.counts, q=0)
beta_div_neutral_treat<-hillpair(data=treat.counts, q=1)
beta_div_phylo_treat<-hillpair(data=treat.counts, q=1, tree=genome_tree)
beta_div_func_treat<-hillpair(data=treat.counts, q=1, dist=dist)
4.5.2.2.1 Richness
Df SumOfSqs R2 F Pr(>F)
time_point 1 2.033529 0.0957379 6.818959 0.001
species 1 2.327685 0.1095866 7.805342 0.001
individual 26 9.722175 0.4577167 1.253885 0.002
Residual 24 7.157208 0.3369589 NA NA
Total 52 21.240598 1.0000000 NA NA
4.5.2.2.2 Neutral
Df SumOfSqs R2 F Pr(>F)
time_point 1 2.168854 0.1067496 9.686182 0.001
species 1 3.090152 0.1520953 13.800733 0.001
individual 26 9.684304 0.4766554 1.663479 0.001
Residual 24 5.373891 0.2644996 NA NA
Total 52 20.317201 1.0000000 NA NA
4.5.2.2.3 Phylogenetic
Df SumOfSqs R2 F Pr(>F)
time_point 1 2.0671823 0.2208963 21.358337 0.001
species 1 0.8731487 0.0933035 9.021460 0.001
individual 26 4.0949687 0.4375828 1.627293 0.008
Residual 24 2.3228577 0.2482174 NA NA
Total 52 9.3581574 1.0000000 NA NA
4.5.2.2.4 Functional
Df SumOfSqs R2 F Pr(>F)
time_point 1 2.1253769 0.4009905 40.2325710 0.001
species 1 0.0070451 0.0013292 0.1333606 0.770
individual 26 1.9000410 0.3584768 1.3833480 0.181
Residual 24 1.2678545 0.2392035 NA NA
Total 52 5.3003175 1.0000000 NA NA
beta_richness_nmds_treat <- beta_div_richness_treat$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(treat_nmds, by = c("sample" = "Tube_code"))

beta_neutral_nmds_treat <- beta_div_neutral_treat$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(treat_nmds, by = c("sample" = "Tube_code"))

beta_phylo_nmds_treat <- beta_div_phylo_treat$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(treat_nmds, by = join_by(sample == Tube_code))

beta_func_nmds_treat <- beta_div_func_treat$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(treat_nmds, by = join_by(sample == Tube_code))

4.5.3 3. Do the FMT work?

4.5.3.1 Comparison between FMT2 vs Post-FMT2

transplant3 <- meta  %>%
  filter(time_point == "4_Transplant2" | time_point == "6_Post-FMT2")

transplant3.counts <- temp_genome_counts[,which(colnames(temp_genome_counts) %in% rownames(transplant3))]
identical(sort(colnames(transplant3.counts)),sort(as.character(rownames(transplant3))))

transplant3_nmds <- sample_metadata %>%
  filter(time_point == "4_Transplant2" | time_point == "6_Post-FMT2")

4.5.3.2 Number of samples used

[1] 53
beta_div_richness_transplant3<-hillpair(data=transplant3.counts, q=0)
beta_div_neutral_transplant3<-hillpair(data=transplant3.counts, q=1)
beta_div_phylo_transplant3<-hillpair(data=transplant3.counts, q=1, tree=genome_tree)
beta_div_func_transplant3<-hillpair(data=transplant3.counts, q=1, dist=dist)
4.5.3.2.1 Richness
Df SumOfSqs R2 F Pr(>F)
species 1 0.9809525 0.0778442 5.775712 0.001
time_point 1 0.7092107 0.0562799 4.175734 0.001
type 1 1.4479860 0.1149060 8.525540 0.001
individual 25 6.9157236 0.5488022 1.628753 0.001
Residual 15 2.5476146 0.2021678 NA NA
Total 43 12.6014875 1.0000000 NA NA
4.5.3.2.2 Neutral
Df SumOfSqs R2 F Pr(>F)
species 1 1.1101711 0.0895457 7.879642 0.001
time_point 1 0.7363935 0.0593971 5.226687 0.001
type 1 1.9027421 0.1534740 13.505059 0.001
individual 25 6.5351386 0.5271203 1.855374 0.001
Residual 15 2.1133659 0.1704628 NA NA
Total 43 12.3978112 1.0000000 NA NA
4.5.3.2.3 Phylogenetic
Df SumOfSqs R2 F Pr(>F)
species 1 0.1341431 0.1005265 6.921553 0.001
time_point 1 0.1159136 0.0868654 5.980943 0.001
type 1 0.1423573 0.1066822 7.345390 0.001
individual 25 0.6512838 0.4880704 1.344205 0.079
Residual 15 0.2907075 0.2178554 NA NA
Total 43 1.3344053 1.0000000 NA NA
4.5.3.2.4 Functional
Df SumOfSqs R2 F Pr(>F)
species 1 -0.0031096 -0.0030953 -0.1177247 0.818
time_point 1 -0.0045864 -0.0045653 -0.1736359 0.827
type 1 0.0775554 0.0771987 2.9361453 0.143
individual 25 0.5385511 0.5360739 0.8155530 0.631
Residual 15 0.3962105 0.3943880 NA NA
Total 43 1.0046211 1.0000000 NA NA
beta_richness_nmds_transplant3 <- beta_div_richness_transplant3$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(transplant3_nmds, by = join_by(sample == Tube_code))

beta_neutral_nmds_transplant3 <- beta_div_neutral_transplant3$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(transplant3_nmds, by = join_by(sample == Tube_code))

beta_phylo_nmds_transplant3 <- beta_div_phylo_transplant3$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(transplant3_nmds, by = join_by(sample == Tube_code))

beta_func_nmds_transplant3 <- beta_div_func_transplant3$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(transplant3_nmds, by = join_by(sample == Tube_code))
p0<-beta_richness_nmds_transplant3 %>%
            group_by(individual) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type, shape=time_point)) +
        scale_color_manual(values=c("#76b183","#d57d2c", "#4477AA")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y = "NMDS2", x="NMDS1 \n Richness beta diversity") +
                theme_classic() +
                theme(legend.position="none")

p1<-beta_neutral_nmds_transplant3 %>%
            group_by(individual) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type, shape=time_point)) +
        scale_color_manual(values=c("#76b183","#d57d2c", "#4477AA")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y = "NMDS2", x="NMDS1 \n Neutral beta diversity") +
                theme_classic() +
                theme(legend.position="none")
  
p2<-beta_phylo_nmds_transplant3 %>%
            group_by(individual) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type, shape=time_point)) +
                scale_color_manual(values=c("#76b183","#d57d2c", "#4477AA")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y= element_blank (), x="NMDS1 \n Phylogenetic beta diversity") +
                theme_classic() +
                theme(legend.position="none")

p3<-beta_func_nmds_transplant3 %>%
            group_by(individual) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type, shape=time_point)) +
                scale_color_manual(values=c("#76b183","#d57d2c", "#4477AA")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y= element_blank (), x="NMDS1 \n Functional beta diversity") +
                theme_classic()+
                theme(legend.position="none")

4.5.4 4. Are there differences between the control and the treatment group?

4.5.4.1 After 1 week –> Post-FMT1

post1 <- meta %>%
  filter(time_point == "5_Post-FMT1")

post1.counts <- temp_genome_counts[,which(colnames(temp_genome_counts) %in% rownames(post1))]
identical(sort(colnames(post1.counts)),sort(as.character(rownames(post1))))

post1_nmds <- sample_metadata %>%
  filter(time_point == "5_Post-FMT1")

4.5.4.2 Number of samples used

[1] 27
beta_div_richness_post1<-hillpair(data=post1.counts, q=0)
beta_div_neutral_post1<-hillpair(data=post1.counts, q=1)
beta_div_phylo_post1<-hillpair(data=post1.counts, q=1, tree=genome_tree)
beta_div_func_post1<-hillpair(data=post1.counts, q=1, dist=dist)
4.5.4.2.1 Richness
Df SumOfSqs R2 F Pr(>F)
species 1 0.6488906 0.0757032 2.115638 0.002
type 1 0.5615418 0.0655126 1.830847 0.006
Residual 24 7.3610760 0.8587842 NA NA
Total 26 8.5715084 1.0000000 NA NA
4.5.4.2.2 Neutral
Df SumOfSqs R2 F Pr(>F)
species 1 0.7984743 0.104299 3.065171 0.001
type 1 0.6051778 0.079050 2.323148 0.004
Residual 24 6.2519772 0.816651 NA NA
Total 26 7.6556293 1.000000 NA NA
4.5.4.2.3 Phylogenetic
Df SumOfSqs R2 F Pr(>F)
species 1 0.0700374 0.0635900 1.6594454 0.140
type 1 0.0184254 0.0167292 0.4365646 0.781
Residual 24 1.0129278 0.9196808 NA NA
Total 26 1.1013906 1.0000000 NA NA
4.5.4.2.4 Functional
Df SumOfSqs R2 F Pr(>F)
species 1 -0.0048931 -0.0058792 -0.1628739 0.837
type 1 0.1161466 0.1395538 3.8660898 0.112
Residual 24 0.7210172 0.8663254 NA NA
Total 26 0.8322707 1.0000000 NA NA
beta_richness_nmds_post1 <- beta_div_richness_post1$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(post1_nmds, by = join_by(sample == Tube_code))

beta_neutral_nmds_post1 <- beta_div_neutral_post1$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(post1_nmds, by = join_by(sample == Tube_code))

beta_phylogenetic_nmds_post1 <- beta_div_phylo_post1$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(post1_nmds, by = join_by(sample == Tube_code))

beta_functional_nmds_post1 <- beta_div_func_post1$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(post1_nmds, by = join_by(sample == Tube_code))
p0<-beta_richness_nmds_post1 %>%
            group_by(type) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type, shape=time_point)) +
        scale_color_manual(values=c("#76b183","#d57d2c", "#4477AA")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y = "NMDS2", x="NMDS1 \n Richness beta diversity") +
                theme_classic() +
                theme(legend.position="none")

p1<-beta_neutral_nmds_post1 %>%
            group_by(type) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type, shape=time_point)) +
        scale_color_manual(values=c("#76b183","#d57d2c", "#4477AA")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y = "NMDS2", x="NMDS1 \n Neutral beta diversity") +
                theme_classic() +
                theme(legend.position="none")
  
p2<-beta_phylogenetic_nmds_post1 %>%
            group_by(type) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type, shape=time_point)) +
                scale_color_manual(values=c("#76b183","#d57d2c", "#4477AA")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y= element_blank (), x="NMDS1 \n Phylogenetic beta diversity") +
                theme_classic() +
                theme(legend.position="none")

p3<-beta_functional_nmds_post1 %>%
            group_by(type) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type, shape=time_point)) +
                scale_color_manual(values=c("#76b183","#d57d2c", "#4477AA")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y= element_blank (), x="NMDS1 \n Functional beta diversity") +
                theme_classic()+
                theme(legend.position="none")
ggarrange(p0, p1, p2, p3, ncol=2, nrow=2, common.legend = TRUE, legend="right")

4.5.4.3 After 2 weeks –>Post-FMT2

post2 <- meta %>%
  filter(time_point == "6_Post-FMT2")

post2.counts <- temp_genome_counts[,which(colnames(temp_genome_counts) %in% rownames(post2))]
identical(sort(colnames(post2.counts)),sort(as.character(rownames(post2))))

post2_nmds <- sample_metadata %>%
  filter(time_point == "6_Post-FMT2")

4.5.4.4 Number of samples used

[1] 28
beta_div_richness_post2<-hillpair(data=post2.counts, q=0)
beta_div_neutral_post2<-hillpair(data=post2.counts, q=1)
beta_div_phylo_post2<-hillpair(data=post2.counts, q=1, tree=genome_tree)
beta_div_func_post2<-hillpair(data=post2.counts, q=1, dist=dist)
4.5.4.4.1 Richness
Df SumOfSqs R2 F Pr(>F)
species 1 0.8966107 0.1132233 3.515588 0.001
type 1 0.6463814 0.0816245 2.534445 0.002
Residual 25 6.3759661 0.8051521 NA NA
Total 27 7.9189582 1.0000000 NA NA
4.5.4.4.2 Neutral
Df SumOfSqs R2 F Pr(>F)
species 1 0.9398968 0.1224370 4.112304 0.001
type 1 1.0227481 0.1332297 4.474800 0.001
Residual 25 5.7139316 0.7443333 NA NA
Total 27 7.6765766 1.0000000 NA NA
4.5.4.4.3 Phylogenetic
Df SumOfSqs R2 F Pr(>F)
species 1 0.1200820 0.1454923 4.647187 0.001
type 1 0.0592745 0.0718175 2.293931 0.027
Residual 25 0.6459931 0.7826902 NA NA
Total 27 0.8253496 1.0000000 NA NA
4.5.4.4.4 Functional
Df SumOfSqs R2 F Pr(>F)
species 1 0.0116738 0.0168395 0.4229082 0.551
type 1 -0.0085273 -0.0123007 -0.3089208 0.909
Residual 25 0.6900904 0.9954612 NA NA
Total 27 0.6932368 1.0000000 NA NA
beta_richness_nmds_post2 <- beta_div_richness_post2$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(post2_nmds, by = join_by(sample == Tube_code))

beta_neutral_nmds_post2 <- beta_div_neutral_post2$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(post2_nmds, by = join_by(sample == Tube_code))

beta_phylogenetic_nmds_post2 <- beta_div_phylo_post2$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(post2_nmds, by = join_by(sample == Tube_code))

beta_functional_nmds_post2 <- beta_div_func_post2$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(post2_nmds, by = join_by(sample == Tube_code))
p0<-beta_richness_nmds_post2 %>%
            group_by(type) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type, shape=time_point)) +
        scale_color_manual(values=c("#76b183","#d57d2c", "#4477AA")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y = "NMDS2", x="NMDS1 \n Richness beta diversity") +
                theme_classic() +
                theme(legend.position="top")

p1<-beta_neutral_nmds_post2 %>%
            group_by(type) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type, shape=time_point)) +
        scale_color_manual(values=c("#76b183","#d57d2c", "#4477AA")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y = "NMDS2", x="NMDS1 \n Neutral beta diversity") +
                theme_classic() +
                theme(legend.position="none")
  
p2<-beta_phylogenetic_nmds_post2 %>%
            group_by(type) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type, shape=time_point)) +
                scale_color_manual(values=c("#76b183","#d57d2c", "#4477AA")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y= element_blank (), x="NMDS1 \n Phylogenetic beta diversity") +
                theme_classic() +
                theme(legend.position="none")

p3<-beta_functional_nmds_post2 %>%
            group_by(type) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type, shape=time_point)) +
                scale_color_manual(values=c("#76b183","#d57d2c", "#4477AA")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y= element_blank (), x="NMDS1 \n Functional beta diversity") +
                theme_classic()+
                theme(legend.position="none")
ggarrange(p0, p1, p2, p3, ncol=2, nrow=2, common.legend = TRUE, legend="right")